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We derive ground state wave functions of superconducting instabilities on the lioneycomb lattice 
induced by nearest-neighbor attractive interactions. They reflect the Dirac nature of electrons in the 
low-energy limit. For the order parameter that is the same irrespective of the direction to any of the 
nearest neighbors we find weak pairing (slowly decaying) behavior in the orbital part of the Cooper 
pair with no angular dependence. At the neutrality point, in the spin-singlet case, we recover a 
, strong pairing behavior. We also derive ground state wave functions for the superconductivity on 

the bilayer honeycomb lattice, with strong interlayer coupling, induced by attractive interactions 
^ ' between sites that participate in a low-energy description. Without these interactions, free electrons 

, are described by a Dirac equation with a quadratic dispersion. This unusual feature, similarly to 

■ ^He - B phase, leads to the description with two kinds of Cooper pairs, with + ipy and Px — ipy 

pairing, in the presence of the attractive interactions. We discuss the edge modes of such a spin- 
singlet superconductor and find that it represents a trivial topological superconductor. 



^ \ I. INTRODUCTION: SUPERCONDUCTIVITY ON HONEYCOMB LATTICE 

o 



X 



The advent of graphene^ opened a door for exploration of new phenomena in two-dimensional Dirac-like condensed 
matter systems. One of the intriguing questions is of superconducting correlations of electrons on the honeycomb 



^ . lattice system. Superconductivity has been induced in short graphene samples through proximity effect with super- 
C/3 ' conducting contacts^. This indicates that Cooper pairs can propagate coherently in graphene. In principle supercon- 
ductivity on the graphene honeycomb lattice can be induced by short-range attractive interactions and explorations 
^ • of allowed possibilities were given in Refsi^"— . Among the most interesting is the so-called p + ip superconducting 
instability introduced in Refj2,. It would be supported by the most natural nearest-neighbor attractive interaction and 



have distinct features of the Dirac electrons. Later it was showed^, by a restricted (low-energy) analysis, that this 
^ I state may be less energetically favorable with respect to Kekule-like order parameter arrangements. Nevertheless, the 
Q . p + ip instability seems, though an exotic state, a very attractive possibility because of its underlying symmetry of the 
O ' order parameter, the same as for Pfaffian quantum Hall stated or p + ip spinless superconductor—. The later systems 
, support non-Abelian statistics, which is at the heart of the idea of the topological computing^. There is an important 
I • difference between these states and the proposed graphene state. The superconducting instability in graphene does 
^ \ not break time-reversal symmetry and those systems do. Due to the valley degeneracy we effectively have two {p±ip) 
\^ • order parameters and that requires additional understanding of intertwined correlations and underlying symmetries. 
' One way, just as in the Pfaffian stated, is to look for the ground state wave function and recognize the structures and 
. symmetries. 

In this paper, in the first part, we will find the effective (long-distance) expression for the ground state wave 
. function of the p + ip spin-singlet instability described in Ref. y and display pertinent symmetries in this case. Also 
' a spinless case will be discussed. We will use the BCS mean-field formalism. In the following section we will set up 
[ the BCS formalism, solve the Bogoliubov - de Gennes (BdG) equations and find the expression for the ground state 
I • wave functions. The last section of the first part is devoted to conclusions. The second part of the paper is devoted 
J> I to the p ± ip superconductivity on the bilayer honeycomb lattice. We refer reader to this part of the paper for an 
introduction. 



II. SUPERCONDUCTIVITY ON HONEYCOMB LATTICE AND ITS GROUND STATES 



The Hamiltonian for free electrons on the honeycomb lattice is 

Ho^~tY^ {al„bj^a + h.c.)-fi'Yfii, (1) 

where t is the hopping energy between nearest neighbor C (carbon) atoms, ai^aiala) on-site annihilation 

(creation) operator for electrons in the sublattice A with spin a ^tiii ^-nd bi rj{bl^) for sublattice B, hi is the on-site 
number operator, and /x is the graphene chemical potential. We use units such that h = I. Diagonalization of Eq.([T]) 
leads to a spectrum given by: = ±i|5'(fc)|, where k is the two-dimensional momentum, and S{k) = ^^exp{ifc(5} 
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with S's defined as Si = a(0,l/-\/3), 62 = a/2(l, — 1/-\/3), and S3 = a/2(— 1, — 1/-\/3), and a = y/3 ttcc, dec is the 
distance between C atoms and a is the next to nearest neighbor distance. At the corners of the hexagonal Brillouin 
zone, K± = (27r)/a(±2/3, 0), we have S{K± + fc) « ^a^/3/2{kx T iky), and the band has the shape of a Dirac cone: 
e{K± + k) = ±VF\k\, where vp = {^/3at)/2 is the Fermi-Dirac velocity. 

For the sake of simplicity we will consider only nearest-neighbor attractive interactions among electrons. The on-site 
repulsive interactions can be introduced and will not change our conclusions. Therefore the complete Hamiltonian 
will include nearest-neighbor interactions as follows, 

where g < 0. We will assume the spin-singlet pairing among nearest-neighbors and apply the BCS ansatz with 
Aij = {tti^^bj^^ — ai_-fbj^^), the superconducting order parameter. Furthermore we assume one and the same Ay = A 
for all nearest neighbors, which due to global gauge (t/(l)) transformations on a's and b's can be chosen real and 
positivo^i. The interaction part, Hj, becomes 

Hbcs = {sE^Kt^L - 4,At^ + - 3<?|A|2. (3) 

The order parameter in the momentum space is 

Ag = ^ A cxp{ik{i~ J)} = A ^ cxp{ik5} = AS{k) (4) 

Therefore near K points A^^^^t ^ ^(kxTiky), which then describes two p- wave like superconducting order parameters 
in a low effective description. The complete BCS Hamiltonian can be now cast in the following form in the momentum 
space. 



where 



k 

k 



= (at ,6t a ,-,,6 ,-, ) (6) 
with defined a^^ = a^cr expjifc i\ and h^^ — hi„ expjifc i}, and, with 5 A = A for short. 



-tS{k) AS'(fc) 

-iS'*(fc) -[I AS{-k) 

AS*{-k) ^ tS(k) 

AS*{k) tS*{k) fi 



We look for the solution in the form of a diagonalized Bogoliubov BCS Hamiltonian, 



Hbcs= E -|.4.."£,.+ E -LC%. + ^o, (7) 

fc,7=± fc,7=± 

where aj: ^ and 7 = ± are new quasiparticles at momentum k. For the dispersions we have: 



wS = 7a;2 and cj5 = , (8) 

fc,7 k fc.7 ' A: ' ^ ' 

where 7 = ±. We define a general solution a as 

Next we have to solve the Bogoliubov - de Gennes (BdG) equations, which follow from the following condition, 

[af:,HBcs]^Eaf:. (10) 
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From this matrix eigenvalue problem we obtain energies of the Bogoliubov quasiparticles, 



En 



±y/ivF\Sik)\+pfj,r + \ASik)\^ 



(11) 



where ± stands for the particle and hole branches respectively for two kinds of excitations p = —1(a) and p = +l(/3). 
For fj, — the system is gapless and we need a coupling g larger than a critical value for the superconducting instability 
to exist^. This can be found considering in the BCS formalism the consistency or gap equation. 

For each valley we have to solve the Bogoliubov problem using the expansion S{K± + fc) « ^a\/3/2{kx T iky). Near 
we need to diagonalize the following matrix, Mi, that comes out of Ea. ([TO)) : 





Vpk 





sk 


Vpk* 


-M 


sk* 








sk 




—Vpk 


sk* 





—Vpk* 


M - 



where s = s* = — Aa-\/3/2 > 0. Its eigenvectors (after normalization) enter the following expressions for Bogoliubov 
quasiparticles: 
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m (m - vpkm^a+^ + b+^) + s\k\U^al^ + b^^)}, (12) 



and 



and quasiholes: 



l^Ep[Ep-{fi + vp\k\) 



:{[Efi- {fj. + vpk)]{^^j —a+-^ - b+-f) - s\k 



(13) 



2^E^[E^ + {ti-vp 



{-[Ea + {}i - vpk)]{\j -^a+T + ^+t; 



^|fc|(\/f7«-i + ^-i)}' 



(14) 



and 



2^Ei3[Ep + ifi + vp 



{-[Ep + {p + vpk)]{^ 



.|fc|(. 



for the Bogoliubov solution near point K^, where we denoted 



K±±k,cr 



a±a and bf% 



K±±kM 



k* 
= b±„. 



(15) 



The natural eigenstates of chirality appeared in our expressions. For example {J ■pra+^ + 6+^) represents spinor: 



(16) 



which is the eigenstate of the chirality operator defined with a ~ {a^, (Jy) Pauli matrices, i.e. the pseudospin (due 
to two sublattices) is along the momentum vector. The state (y^^O-i + ^-^) represents the same spinor because of 

the interchanged roles of sublattices in the A'_ point. To see this in more details we would like to remind the reader 
that instead of the Dirac free electron representation by the spinor 



'^k ^ K++k,a' K++k,a K-+k,a' K-+k,a" 



(17) 



and the chirality operator is defined as 



3k Q 

lfcl 



\k\ 



(18) 
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in the BdG formalism wc work with 



^ X+ + fcT K++kt K_-kl^ K^-ki' 



[a 



(19) 



Note the reversed order of sublattices and the change of the sign of the momentum k near point in the BdG 
formahsm with respect to the free one. Thus the lower 2x2 matrix on the diagonal of the Hamiltonian matrix in the 
free Dirac case can be read off from: 



at . 

K^+k,a K-+k,a 



—fi —vpk* 



^K-+k,a 



(20) 



i.e. it is equal to —vpka — ji. Note that if we change the sign of k vector in Eq. (|20l) i.e. k — >■ —k the off-diagonal 
elements in the matrix will change the sign, so that in this basis in the free representation the chirality operator will not 



have minus sign in the lower right entry of the matrix representation in Ea. (jl8|) . Therefore [y 



-i' 



-&_^) represents 



the same spinor (up to a phase factor) as in Ea. (|16p and the same chirality eigenstate (with positive eigenvalue) as 
we pointed out earlier. Nevertheless in the Bogoliubov representation wc still have 



-vpk 



-vpk* 



(21) 



i.e. the matrix is —vpka + /i, and the representation of the chirality operator stays the same as in Eg. psp . We will 



use this fact later. On the other hand the combinations in Eqs. and (|57|: (w p-a+^ 

have the pseudospin vector in the opposite direction of the momentum vector k. 
It is thus natural to introduce the following notation: 



5+t = C+tu., 



- J 



-b+t) and (y'^a_^-fe_;) 

(22) 
(23) 
(24) 
(25) 



where v and w denote the chirality i.e. whether the pseudospin vector is along or in the opposite direction with 
respect to the k vector, respectively. We have to note that these electron operators are defined up to a phase factor, 
most importantly phase. This degree of freedom should not influence the physics, but we chose the definitions 
so that later the symmetry under exchange of particles in the ground state wave function is transparent. 

The a and /3 sectors are obviously decoupled in the Bogoliubov description and we can concentrate and closely 
examine the a sector first. Furthermore we do not have to consider K_ point separately as the symmetry considerations 
tell us that the BdG equations around this point will induce the coupling or states of an electron around A'+ point 
with \. projection of spin and those around K_ point with 'f projection of spin. 

Thus it suffices to consider a sector first (with c+^t, and c_^i,) and then use the symmetry arguments, more precisely 
antisymmetry under real spin exchange to recover the whole ground state wave function. Wc can rewrite a's in the 
following form. 



ak- u'^c+^v + «fccl^„. 



(26) 
(27) 



We should demand ak.+ \G) = and a^. _\G) = 0, for any k, if \G) is to represent the ground state vector. That 
implies that in the a sector of point we have the following contribution to the ground state. 



n 



(28) 
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where |0) denotes the vacuum. This state is annihilated with both, ak,+ and a\. _. The symmetry arguments demand 

that we should get a similar expression considering BdG equations at point. If we denote by ga{k) = — -f, the 
ground state vector in the a sector should look like: 

k 

= nil + 9.ik)[ci^j_,^ + ci^j^j + ^[i,„ct_^, + 

k 

= exp{^g„(fc)[c;^,cl^„+cl^„c^+^j}|0) (29) 

Now we can identify ga{k) to represent a Fourier transform of the wave function of a Cooper pair of electrons, which 
is a spin-singlet with respect to spin degree of freedom and a triplet state (symmetric under exchange) with respect 
to valley {K±) degree of freedom. If we defined differently our electron operators there would be possibility for ga{k) 

to acquire the phase factor ^J^, which would make the identification of the antisymmetry under exchange harder. 

Taking into account the /3 sector (with the chirality in the opposite direction of the momentum: w) the complete 
ground state vector is 



exp{^ 



g„(fc)[c;^,ct_^, +ct_^,„ct^,,J +^.g,(fc)[ct^^^ct_^^ +ct_^^„ct^,^]}|0), (30) 

k 



where 



s\k\ s\k\ 
ga{k) = -— r-jr and g/}{k) = -— — ■ — -. (31) 

Ea - [fl - VF\k\) Ea - [fl + VF\k\) 

Using the long-distance (low-momentum) expansions for Ea and Ep, for finite /x, 

Eo.ip)-I^TVF\k\ + ^jj-, (32) 
we find the long-distance behavior of the pair wave function to be 

lim ga{r) ^ \im .g/j(r) - j^r, (33) 

|r|->oo |r|-i.oo \r\ 

i.e. we have a case for a weak pairing^. As emphasized in Ref. the term weak pairing does not mean also weak 
coupling, it stands for a phase with an unusual large spread of the Cooper pairs. On the other hand for /i = we 
have that ga{k) and gp^k) are two constants and the Cooper pairs are localized on a short scale '--^ a in the graphene 
system at the neutrality point. Thus for fi ~ wc have a case for a strong pairing. 

The ground state vector (wave function) in Eg. pop displays two kinds of Cooper pairs, each antisymmetric under 
combined exchange of (a) orbital, (b) valley {K±), and (c) spin (t,i) degree of freedom. Two kinds of Cooper pairs 
stem from the chirality (sublattice) degree of freedom intimately connected with the Dirac-nature of the electron with 
both, particles and holes. They both, particles (with positive chirality v at Kj^) and holes (with negative chirality w 
at K^), constitute Cooper pairs, which are symmetric under v w,vf —vp transformation. 

In the long distance limit we recover the form of the wave function of ordinary s-wave superconductor as given in 
Ref. |iq[ though with more, two-component, degrees of freedom. The Cooper pair wave function is antisymmetric 
under spin exchange and symmetric under exchange of valley {K±), sublattice (u, w), and orbital degrees of freedom. 

Next we will discuss the spin-triplet case, more precisely we will assume that the system is spin-polarized and not 
consider spin in the following. Therefore fermions are spinless just like in the Pfaffian case, but they live on the 
honeycomb lattice. We will assume {aibj) = A. In this case the Bogoliubov problem in Eq.(IS|) for the spin-singlet 
pairing transforms into a similar one with ^ = and ^ = bj:, and the matrix becomes as follows 



-/i -tS{k) AS{k) 

-tS*{k) -AS'(-fc) 

-A5*(-fc) fi tS{k) 

AS*{k) tS*{k) ^ 
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Around the Kj^ point we have 



—/.J vpk* sk* 
vpk — /i —sk 
—sk* fi —vpk* 
sk —Vpk fi 

where s = —Aa^ > as before. The problem around the /v_ point is a copy of the problem around the Kj^- point. 
Now the Mj: matrix around point cannot be cast, as in the spin-singlet case, in the following form, 



vpdk — /1/2 sak 

sak —Vpdk + i^il2 



where I2 is the 2x2 identity matrix, which commutes with the chirality matrix fEallSp. around /^_|_ point can 
be compactly written as 



Vpdk — /i/2 sik X d 
—sikxd —vpdk + fil2 

and it does not commute with the chirality operator. The eigenstates of the Bogoliubov problem do not have to be the 



eigenstates of chirality. We find the following eigenvalues Ep — ±y fi"^ + \k\'^s'^ + \k\'^Vp + p 2\J ii^v'^\k\'^ + s2?j|,|fcp, 

where p = +1(0;) and p — — are two branches as before. The associated eigenvectors can be written as sums of 
fermionic particle eigenstates of chirality only in the low-momentum limit and we list those connected with positive 
eigenvalues. 



2(1 



(34) 



and 



2(1 



|fe|2£2 , 



2fi 



(35) 



and negative eigenvalues. 



2^ 



2(1 + ^) "^'^ 



(36) 



and 



fe,- 



2(1 



4^^ 
\k\^s 



(37) 



Similarly as before we can define 



k* 



- J 



(38) 
(39) 
(40) 
(41) 
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and the ground state vector can be cast in the fohowing form, 

^^P{E + 4.»c^-,.)}|0)- (42) 

In this case each Cooper pair is antisymmetric under exchange of K± points i.e. vahey degree of freedom and 
symmetric under exchange of sublattices i.e. chirahty (w o w). Depending on our definitions for c's two degrees of 
freedom can exchange the symmetry properties. We find again the weak pairing (~ i) behavior in the orbital part. 



III. CONCLUSIONS: SUPERCONDUCTIVITY ON HONEYCOMB LATTICE 

We derived the ground state wave functions for the superconductivity on the honeycomb lattice induced by nearest- 
neighbor attractive interactions and with order parameter independent of the direction to any of the nearest neighbors. 
Although the order parameter in momentum space has the p±ip form in a low effective description the Cooper pair 
wave function behaves as s-wave (with no angular dependence) and decays as ~ -. Other (discrete) degrees of 
freedom combine to make the Cooper pair antisymmetric under exchange. At the point of the transition, yu = 0, in 
the spin-singlet case a strong pairing (of the order of lattice spacing) occurs. 

IV. INTRODUCTION: SUPERCONDUCTIVITY ON BILAYER HONEYCOMB LATTICE 

Topological superconductors in a strict sense or what we also call non-trivial topological superconductors have odd 
number of Majorana modes moving in each direction on the edge of such a superconductor—. In the case of trivial 
topological superconductors we have even number of Majorana modes i.e. by combining them in pairs we can talk 
about Dirac fermions on their edge, s-wave superconductor in two dimensions is always topological in the sense that 
it has a gap in its bulk and non-trivial degeneracy of the ground state on the torus (equal to four)^. We have to 
use one Bose field (one Dirac fermion) to describe the edge of such a systenJ^. On the other hand if we combine two 
p-wave superconductors, with + ipy and Px — ipy orbital symmetry, and each of the two corresponds to one let's 
say definite projection of spin we have the case for a non-trivial topological superconductor. In that case on the edge 
live two Majorana modes that are moving in opposite directions and each is associated with different projection of 
electron spin. This represents a "helical" edge where we have a pair of edge Majorana modes (moving in opposite 
directions) that are connected with a time-reversal operation. This is the simplest topological superconductor we can 
imagine in two dimensions and has yet to be realized and detected in experiments. In three dimensions a realization 
of topological superconductor is He^ B - phase^^. 

On the other hand the honeycomb lattice, nowadays very much connected with the research on graphene, is a 
playing ground for various, among others topological, phases. The first topological insulator was introduced on the 
honeycomb lattice with a special interactioni^. While considering possibilities for superconducting instabilities on the 
honeycomb lattice and graphene in Ref. [3, a phase was proposed with two p-wave order parameters (each near two 
effective descriptions in fc-space i.e. two valleys). Though one might expect that, while considering triplet pairing i.e. 
if we suppress spin, this would lead to a non-trivial topological superconductor with a pair of Majorana modes, this 
is not the case as we demonstrated in the first part of the paper (Ref. IT^ . We found that the ground state wave 
function is antisymmetric with respect to the valley degree of freedom and that the orbital part of the Cooper pair 
is with no angular dependence i.e. a s-wave. As we emphasized earlier in the case of s-wave we expect one Dirac 
fermion per degree of freedom on the edge i.e. no non-trivial behavior. 

In this paper we will derive ground state wave functions for some superconducting instabilities that may emerge 
due attractive interactions on a two layer system in which each layer represents a honeycomb lattice. The two lattices 
are stacked as in the bilayer graphene i.e. in the way of Bernal stacking. As in the bilayer graphene, we expect that 
the pseudospin vector connected with sublattice degrees of freedom will not follow the momentum vector in a parallel 
or antiparallel fashion, like on ordinary honeycomb lattice or graphene, but rotate for a whole angle as a result of a 
rotation of the momentum vector for a half an angle. This feature of free electrons on the bilayer honeycomb lattice 
will refiect in the description of Cooper pairs when attractive interactions are introduced. The explicit ground state 
wave functions and Cooper pair structure will help us to see more closely the nature of pairing in this system. We 
will find the p-wave angular dependence in the orbital part. 
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B1 




81 



B2 



FIG. 1: A view of Bernal stacked honeycomb lattices 1 and 2 with corresponding sublattice sites Al and Bl, and A2 and B2. 

In the following section, we will formulate Bogoliubov - de Gennes (BdG) equatons for this system. In the next 
section the explicit solutions with corresponding ground state wave functions will be given in the case of (a) spin- 
singlet and (b) spinless (spin-triplet) pairing. Then we will examine whether these systems arc truly gapped in the 
bulk, and, in the spin-singlet case, its edge spectrum. The last section is devoted to discussion and conclusions. 



V. ELECTRONS ON BILAYER HONEYCOMB LATTICE AND BCS INSTABILITY 



The Hamiltonian for free electrons on two honeycomb lattices, which are Bernal stacked, is 

Ho = -tY^ ^ia.l,,Xah,n+la + 0^2M,aKn-S,a + '^'C-) + ^-L ^1 («ln,^«2,S,^ + h.c) - ^l^nn. (43) 

The index i = 1,2 denotes the layer index. In Fig. 1 the relative positions of two triangular sublattices, Ai and 
Bl, for the lattice 1, and A2 and B2; for the lattice 2 are illustrated. In Ea. (P5|) t is the hopping energy between 
nearest neighbor C (carbon) atoms in the case of the bilaycr graphene in each layer, and t±^ is the same energy for 
hopping between the layers. The on-site creation (annihilation) operators, a\ fi ^{cLifi^a), are for the electrons in the 
sublattice Ai of the layer i with spin a =t,i, and 6]^ - ^{bi^n.a) for the electrons in the sublattice Bi, fifi is the on-site 
number operator, and /i is the chemical potential. (5's are defined as 5i ~ a(0,l/A/3), ^2 — a/2(l, — 1/\/3), and 
^3 = a/2(-l, -l/\/3), and a V3 flee, flee is the distance between C atoms and a is the next to nearest neighbor 
distance. 

We use units such that h=l. By introducing Fourier transforms ^ = X)n '^i,",'^ ^-'^Pi*^ ^} ^^'^ ^ika ^ 
'^fibifi^c! cyi'p{ik n\ etc. and diagonalizing the Hamiltonian we find for the spectrum. 



E^ik) = ±{{-\r\ + t^\s{m , (44) 

where S{k) = ^^exp{ifc 5}, and a = 1,2 stand for two kinds of branches. Near K points, the corners of the hexagonal 
Brillouin zone, K± = (27r)/a(±2/3, 0), we have 

S{K± + fc) « Ta\/3/2(fc, T iky), (45) 
and in the limit tj^ ^ t the lower positive and higher negative branch have the folowing dispersion relation, 

E^^±^, (46) 

1 2m* ^ ' 
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where m* = and vp = (A/3ai)/2, the Fermi-Dirac velocity. The effective Hamiltonian near K pointSiiSi is 



Hef(k)=-f 





(fc)2 



and acts on the subspace of (pseudo)spinors 



a{K++k)i 



around point and 



around point . Hgf can be rewritten as 



Hefik) 



1 

^ "1.2 u2\ ^. . , K 



2m 



2m 



-an, 



(47) 



(48) 



(49) 



(50) 



where k — |fc|(cos{0j}, sin{0^:}) and n — (cos{20^}, sin{20^}), and cr's are Pauli matrices. The operator an encodes 
the projection of the pseudospin on direction n. For eigenstates as 



X 



, and X 



k 

-1 



(51) 



the direction n may be interpreted as the direction of the pseudospin vector, with projection (chirahty) equal to +1 in 
the case of x+ , and — 1 in the case of x~ ■ Thus in the case of these eigenstates we see explicitly our previous remark 
that the pseudospin vector rotates for an angle while k vector rotates for half an angle circling the Fermi surface 
around K points. That feature of the solutions of the free problem leads to non-trivial pairing in the orbital part 
of Cooper pairs as we will see later. This is to be contrasted to the behavior in the monolayer, a single honeycomb 
lattice, where the rotation of k vector is strictly followed by the rotation of the pseudospin vector. It is accompanied 
by s-wave pairing, when special (nearest-neighbor) attractive interactions are applied. In that case although two order 
parameters are of, + ipy, and — ipy type we have the trivial (s-wave) behavior in the orbital part as we have 
shown earlier. 

As the reader may have noticed we did not include the direct hopping between the atoms of Bl and B2 sublattice. 
This inclusion is required when we model bilayer graphene^^, but even there for realistic parameters this does not 
influence the physics at high electron momenta or strong magnetic fieldsi^. 

But we will consider nearest-neighbor attractive interactions between electrons on Bl and B2 sublattice. Namely 
these sublattices by themselves make a honeycomb lattice as we can verify by looking at Fig. 1. Due to the strong 
hopping between Al and A2 sublattice the complete low-energy physics is projected onto Bl and B2 sublattice. If the 
interactions are not too strong they can be simply added to this low-energy subspace. The on-site repulsive interactions 
can be introduced and we do not expect that will change our conclusions. Therefore the complete Hamiltonian will 
include nearest-neighbor attractive interactions between electrons on Bl and B2 sublattice as follows, 



(52) 



where g < 0- We will assume the spin-singlet pairing among nearest neighbors and apply the BCS ansatz with 

the superconducting order parameter. Furthermore we assume one and the same = A for all nearest neighbors, 
which due to global gauge (U(l)) transformations on 6i's and &2's can be chosen real and positive. The interaction 
part. Hi, becomes 



Hbcs = {^E ^(^l.- T^l.+u - <r.An^i,) + h-c} - 3,|A| 



ft, 5 



The order parameter in the momentum space is 

Ag = E A exp{ik5} = AS{k). 



(54) 



(55) 
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Therefore near K points A^^^^r '-^ T{kx T iky)^ which then describes two p-wave like superconducting order param- 
eters in a low-energy effective description. Taking into account the complete low-energy reduction the total BCS 
Hamiltonian can now be cast in the following form in the momentum space near K^, q = + /c, 



Hbcs = J2 ^l^'k^?^ 



(56) 



where 



and, with s = s* = —Aga\/3/2 > 0, 



4 = ^U-T^2,-g-i, 







{k'f 

2m* 





sk* 




_ IMl 

2m* 




sk 




2m* 







sk* 






sk 





2m* 


M - 



(57) 



We will omit the discussion concerning Mq' in the neighborhood of A'_ and momenta: q = — k. This entails 
operators which combine f spin with momenta q = K- — k and \. spin with momenta q = + k, and will not provide 
any new information for the structure of the ground state wave function or energy dispersion at small momenta. We 
can simply include these operators at the end in the ground state wave function following symmetry requirements for 
the spin-singlet pairing. 



VI. GROUND STATE WAVE FUNCTIONS OF SUPERCONDUCTING INSTABILITIES 

Wc look for the solution of Eq. ((56|) in the form of a diagonalizcd Bogoliubov BCS Hamiltonian, 

Hbcs - E ^lA..^t.+ E A. + ' (58) 



:,7 fc,7 ''■'T ^ — ' fc,7 fc,7 '^'T 



where ctj:^ and /3jr^, 7 = ± are new quasiparticles at momentum k. For the dispersions we have: 

uj'i = -iuj'i and uji = , (59) 
where 7 = ±. We define a general solution a as 

Next we have to solve the Bogoliubov - dc Gennes (BdG) equations, which follow from the following condition, 

[ttj:, ifscs] = Eaj:. (61) 
We need to diagonalize the following matrix, MZ, that comes out of Eq. (|6ip : 





(kf 

2)ri* 





sk 


2rn* 




sk* 








sk 




2m* 


sk* 





(fc*)^ 

2m* 


A* - 



(62) 



From this matrix eigenvalue problem we obtain energies of the Bogoliubov quasiparticles, 

^ = ±^M^ + \WS^ + (^)'l^l' +?'y4(^)2|fc|4 + 2{^m^ - (fc*)2|fc|4(_f_)2 _ (fc)2|fc|4(_A_)2^ (53) 
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where ± stands for the particle and hole branches respectively for two kinds of excitations p = —1(a) and p = +l(/3). 
The eigenvectors of matri: 
Bogoliubov quasiparticles 



The eigenvectors of matrix Mi (after normalization) enter the following expressions in the long-distance limit for 



and 



and quasiholes: 



and 



= 7f {(|^)(^2+t - h+t) i^bl_^ bl_^)h (65) 
".V = + + + bU)h (66) 

- ^{-i^b^n - ^i+t) - (0)(4-, - bU)}, (67) 

for the Bogoliubov solution near point K^, where we denoted Kj_±k a ^ ^2±(t SiJid K^±k a ^ b\±„. 

It is helpful to introduce the creation and annihilation operators of states of definite chirality or projection of the 
pseudospin along the n vector; we denote by v positive projection, and by w negative projection: 

(68) 
(69) 
(70) 
(71) 

Notice that the states are defined up to a phase factor according to the definitions in Ea. ([5T|) . and that in the states 
around /v_ point the roles of electrons in different layers are interchanged. We may then define their superpositions, 

C±r = (c±« + c±^)/2, (72) 
c±i = (c±^ - c±„)/2. (73) 



C+i, = 


k^b2+t 


+ bi+^, 


C-f — 


k* 


- bi+-t, 


C—u — 


k* 


+ bi-i, 






4 + bi- 



Then we can rewrite Bogoliubov operators as: 



-^{rkc+i + rk* c+r + cL^ + cl^} , (74) 
V2 



I3k,+ = —{~rkc+i + rk*c+r + cLr - c_;}, (75) 



1 

ctk- = —^{rkcL^ + rk*c}_i - c+r - c+(}, (76) 
1 



Pk- = -^{rfccLr - rk*c^_i - c+r + c+i}, (77) 



with r = ^-nd find that the wave function, 5*0, 



*0 = (1 - i^rcU - i-Aicl, + -2kl2^+r^-Al^-M (78) 



rk* rk 

is annihilated by quasiparticle annihilation operators, a^^+^'o = fik,+^o = 0, and by quasihole creation operators. 



a\ _\^o = - ^13 = 0. Therefore the new ground state is 



k k 
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where we explicitly introduced the sector that couples momenta around with spin J, and momenta around K- 
with spin f by enforcing the explicit spin-singlet pairing which we introduced at the beginning. 

From the structure of the ground state wave function for spin-singlet pairing in Eq. (|79| we find that electrons pair 
between and K- point have the same pseudospin and therefore we have two distinct Cooper pairings for two 
orthogonal pseudospin states, which we denoted by r and I. Each pair has a p-wave pairing in the orbital part and, as 
we work with a time-reversal invariant system, two distinct Cooper pairings are accompanied by two distinct, + ipy 
and px — ipy-, symmetries in the orbital part. Each Cooper pair is antisymmetric under spin exchange, valley exchange, 
and exchange in the orbital part and symmetric under sublattice (pseudospin) exchange. Therefore the ground state 
wave function in Eq. ((79)) describes a Cooper-paired collection of fermions -electrons on the bilayer honeycomb lattice 
with unconventional p-wave pairing. 

The non-trivial (non-s-wave) pairing cannot be eliminated with a gauge transformation. To preserve the form of the 
free part of the Hamiltonian any gauge transformation should be fe^^^ — )■ exp{i0+}62+a- ^^'l ^i+cr ^ 6xp{i0-)_}6| 



and b 



exp{i(/)_}fe2_cr b\_^ — > exp{i0_}6|_jj, and if we rewrite in terms of these operators 

= exp{- i:(blA-i blA-t) E 



l + CTJ 



(80) 



we see that by this gauge transformation we cannot eliminate simultaneously the angular dependence in the two types 
of Cooper p-wave pairings. 

Next wc will discuss the spinless case. We will assume that all electron spins are polarized and that (&in&2fi+(5) ~ 
In this case the Bogoliubov problem in Ea. (|56p for the spin-singlet pairing transforms into a similar one with b^j:^ = b^^ 



and b^j:^ = b^^ and the matrix becomes as follows. 



Qiq—K^+k k 





2m* 





-sk* 


{kf 

2 m* 




sk 








sk* 




(k'f 


—sk 





2in* 





The problem around the point is a copy of the problem around K^. We find the following eigenvalues for the 
Eq.dnil), 



^ - ±^A^^ + + (ir)'l^l' +?'\/4(^)2|fc|4 + 2{^m^ + (fc*)2|fc|4(_f_)2 + (fc)2|fc|4(_i_)2^ (gl) 

where ± stands for the particle and hole branches respectively for two kinds of excitations p = — 1(q;) and p = +1(/?). 
The eigenvectors enter the following expressions in the long-distance limit for Bogoliubov quasiparticles : 



and 



and quasiholes: 



and 



a 



Introducing as in the spin-singlet case the following pseudospin operators 

J -At 



cl; = 7f4- ^nd cL^ = b\_ 
x4+ and = b\^ 



(82) 



(83) 



(84) 



(85) 



(86) 
(87) 
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we can rewrite the eigenvectors as follows, 





—^{rkc+i - rk*c+r + 






(88) 




— {-rkc+i - rk*c+r 


+ cL„ 




(89) 




—={rkc' ,,, — rk* c' , 4 


C+r 4 


-c+i}, 


(90) 


Pk- = 


—={rkcLj, + rk*c^_i + 
v2 


C+r - 


c+;}- 


(91) 



Similarly as in the spin-singlet case we can find that the ground state wave function can be expressed as 

A; fe 

= ^-piE ^biA- E ^bi_bUm- (93) 

A; k 

The Cooper pairs are symmetric under valley and sublattice (pseudospin) exchange and antisymmetric under exchange 
in the orbital part. We have two kinds of Cooper pairs with underlying px + ipy and px — ipy symmetry. 

VII. THE NATURE OF PAIRING PHASES 
A. Gaps and possible nodes 

In Ref. [13 the case of spin-singlet pairing in the monolayer was thoroughly discussed. A topological phase structure 
was described with four (due to valley and spin) Dirac edge modes. This is consistent with our previous calculations 
of the ground state wave function that has explicit s-wave dependence ~ for |^| > (chemical potential). Here we 
extended ground state calculations to the spin-singlet and spin-triplet case of the bilayer. It is appropriate to ask the 
question raised in Ref. [13 for spin-singlet and spin-triplet monolayer case: Are these phases truly gapped or there are 
nodes for some fc's in the bulk spectrum ? In the same reference it was found that in the spin-triplet case, as opposed 
to the spin-singlet case, there are nodes in k space at which gap is equal to zero. We will find a similar situation in 
the bilayer case, except that in the case of spin-singlet pairing there is a critical value for chemical potential above 
which we have a truly gapped - topological phase. (The triplet case just as in the monolayer analysis in Ref. [l3 has 
nodes in the bulk spectrum.) 

It is not hard, by repeating the approach of Reffisl. to find expressions for the matrix, Alp that enters BdG 
equations in both cases for general (not low) momentum k. They are 



Mr 

k 



-M -T{S*{k)f ^S*{k) 

-T{S{k)Y -n ±ASik) 

±AS*{k) fi r(S'*(fc))2 

ASik) T{S{k)f n 



(94) 



where T = j— and + and — stand for spin-singlet and spin-triplet case respectively. We look for the zeros of Bogoliubov 



quasiparticle energies, expressed in the low k limit in Eq. (|63|) and Eq. ()8ip . when p = —1. With no low momentum 
limit, from Ea. (|94[) . their expressions are, 

= ^J^l■2 + |5(fc)|2A2 + \S{k)\^v^ - ^4\S{k)\^^i-^T2 + A2T2|5(fc)|4[2|5(fc)|2 ^ 5(fc)2 ^ (S^ik)^], (95) 

where we used shorthand notation taking gA = A, and E~ and correspond to the spin-singlet and spin-triplet 
case respectively. We assume that A is small with respect to fi so that possible nodes can be only near Fermi surface 
defined by = r|5(fc)|2. Equating E~ and £'+ to zero is equivalent to the following condition, 

- T^\S\^f = -I^I^A^ - 2/i2|5|2A2 T 2A'T^\Sf cos{2(/)}, (96) 
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where is the phase of S{k) - complex number in generaL If we assume that we work with fc's near Fermi surface we 
have approximately 

2//(l ± cos{2(/)}) « T4rcos{2(/)}5 - ^A^, (97) 

where S is defined by |sp — -yr+J as a small depature from the Fermi surface value. Therefore we can approximate 
that for possible nodes near Fermi surface in the case of spin-singlet pairing S{k) is imaginary i.e. (j) — ±^ and in the 
case of spin-triplet pairing S{k) is real i.e. = 0, tt. 

It is not hard to find nodes in the spin-triplet case. From the definition, 

S{k) = exp{z^} + exp{zi(fc. - ^)} + exp{-zi(fc. + -^)}, (98) 

we can recognize that possible positions of nodes can be restricted to fc^; = because in that case S{k) takes real 
values; S{k) = 1 -|- 2cos{-^}. Then in our approximation from Eg. ip?)) and definition 15*^ = ^ + 6 we have 

ky = \/3arccos{^^^ ^ ^J" ^ a/3 arccos{ ^ ^^^'^ ^^^^ 

k 

and that with kx = ^ defines a position of a single node. Therefore the existence of this node (and other related by 
symmetry) tell us that this phase is likely to be gapless even in the bulk and can not represent a topological phase. 

In the spin-singlet case, if we rescale the momentum fcy as — ky in Ea.(|M)). the condition that S{k) is purely 
imaginary demands that 



k k 

cos{fc.y} + 2cos{^}cos{^} = 0. (100) 



Then 



2/m5(fc) = sin{fc,J _ 2 cos{^} sin{^} = ""'"^ f\ (101) 

2 2 (2os{-^} 

Expressed differently as 

(1 - cos{3fcy}) = 4|5p(l + cos{ky}), (102) 

this leads to the conclusion that for large enough |5p i.e. chemical potential this equation does not have a solution 
for cosjfcj,}. We find that for jS'p > |(a/T2 — 3) ~ 0.348 no solution exists. Therefore for large enough chemical 
potential we can have a spin-singlet topological phase i.e. a phase with no gapless bulk excitations. 

B. Edge modes 

To further examine the topological nature of the spin-singlet phase we will derive its edge modes. With respect 
to the lattice structure we will consider a particular geometry where the system, defined on a half-plane, has the 
edge at X = 0. We remind the reader that we use the convention in which K± vectors are along x axis. This choice 
of boundary corresponds to so-called armchair boundary condition for which we require that the solutions of BdG 
equations vanish at x = 0. 

We will consider BdG equations in the low k limit around K± points, neglect terms quadratic in k{k*), and employ 
the substitution k^ — > — to get their form in the real space; due to the symmetry of the problem we seek solutions 

in the form cxTp{ikyy} f (x) and keep the ky dependence. The expression for BdG matrix at momentum + A; is 



-^0 sk 

-Ai sk* 

sk ^ 

sk* ^ 



(103) 
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and at momentum — fc, 





sk 
sk 



sk* 

~fi sk 

^l 

^ 



(104) 



and wc immediately see a reduction of the problem to four sets of equations with the substitution k^ 
denote by 

f{x) = {u^,ui,vi,vi)'^, 



If we 

ox 



(105) 



a form of a solution, where we used ^ and ^ to denote in a shorthand notation pseudospin (compare with Eq. (j60p 
where the same symbols were used for real spin), we have four sets of equations of similar form. For example, one set 
at point is: 



d 

- /.m-f — is—v-f + iskyV^ = Eu^ 
d 

IIV+ — is—Uf — iskyUf ~ Evf. 
ox 



(106) 



Just as argued in Ref. in the case of a simple p-wave superconductor, when E — and fc^ = we have a zero 
mode, u-f = iv-f ~ exp{— -^a;} exp{«^}. The difference here is that and carry opposite valley indexes and we 
cannot construct Majorana Bogoliubov quasiparticle with them, but only a Dirac one. For E we find a chiral 
(unidirectional) mode: E = sky,ky > 0. Considering analogous equation for u-f and at point we get an 
additional solution with the same chirality (the same direction of ky), but these two solutions enter the boundary 
condition which requires that vanishes at the boundary. Therefore one solution with fc^ > is 



explikyi/} exp{ x} exp{i — } sinjQx}, v-f 



-iuf, = = 0, 



(107) 



where Q = |A'-|-|. The other, with ky < and E = —sky, we find in analogous way. Therefore, we find, as expected 
from the ground state wave function in Eq. (j79p . two (Dirac) solutions with opposite directions of motion along the 
edge that carry opposite pseudospin. In addition, due to the spin degree of freedom we expect doubling of modes i.e. 
in total two Dirac modes in each direction. 



VIII. DISCUSSION AND CONCLUSIONS: SUPERCONDUCTIVITY ON BILAYER HONEYCOMB 

LATTICE 



We do not have to go into a discussion of topological invariants to recognize that ground states for spin-singlet pair- 
ing, Eq. ([7^ . and for spin-triplet pairing, Ea. (|M|) . may represent ground states of trivial topological superconductors; 
spin and valley degree of freedom induce the doubling of Majorana modes that may follow^ from p-wave pairing in 
the orbital part. Nevertheless these models of superconductors are interesting in their own right due to the presence 
of the p-wave pairings in the orbital part. In the case of spin-singlet pairing we found that there are no gapless bulk 
excitations for large enough chemical potential and this phase can represent a trivial topological superconductor. 

In deriving ground state wave functions, in the spin-singlet and spin-triplet case, we assumed validity of small k 
expansion around fc = 0. In a strict sense this requires that kp is in the same neighborhood where we can approximate 
dispersion relations i.e. either effective hopping (Fermi velocity) is large or chemical potential is small. We also allowed 
the possibility that minima or even nodes (in the spin-triplet case) in the spectra of the superconductors can be away 
from fc = 0, around kp. In the cases of the trivial topological superconductors, on the monolayer and bilayer 
honeycomb lattice, this seems completely justified in the view of their edge spectrum (see alsoii). In the cases of 
gapless spin-triplet superconductors derived ground state wave functions, though we are inclined to associate them 
with topological phases, may describe these critical statesi^. We find that only if we apply bias to the bilayer larger 
than its chemical potential we can not apply the BdG program around fc = in the way we described in the paper. 
(We remind the reader that with the bias the energy minima of the noninteracting problem shift from fc = at A' 
points to nonzero fc's.) 

The attractive interactions that we need for the realization of the paired ground states and phases may well be 
within the reach of future experiments. In the case of a single honeycomb lattice the interactions may be induced by 
chemically doping the graphene via metal coating^ or trapping fcrmionic atoms in a honeycomb optical lattice^^. 
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To conclude, in the second part of this paper we derived ground state wave functions for the superconductivity on 
the bilayer honeycomb lattice (with strong interlayer coupling) induced by attractive interactions between sites that 
participate in a low-energy description. As is well-known, without these interactions, free electrons are described by a 
Dirac equation with a quadratic dispersion. This unusual feature, similarly to '^He - B phase, leads to the description 
with two kinds of Cooper pairs, with + ipy and Px — ipy pairing, in the presence of the attractive interactions. This 
is expressed in Eq .([79|) in the case of the spin-singlet pairing. Due to the spin degree of freedom we find doubling of 
two chiral Dirac modes with opposite pseudospin on the edge of this spin-singlct superconductor - a trivial topological 
superconductor . 
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